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Abstract 

Using a simple hydrophobic/polar protein model, we perform a Monte Carlo study of the ther- 
modynamics and kinetics of binding to a target structure for two closely related sequences, one of 
which has a unique folded state while the other is unstructured. We obtain significant differences 
in their binding behavior. The stable sequence has rigid docking as its preferred binding mode, 
while the unstructured chain tends to first attach to the target and then fold. The free-energy 
profiles associated with these two binding modes are compared. 
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I. INTRODUCTION 



Protein structures are often viewed as well-defined static entities. For many proteins, 
this simplified static picture is accurate enough to provide valuable information about the 
function. 1 However, there are proteins that are wholly unstructured or only partially struc- 
tured and yet functional. In fact, recent studies suggest that such proteins are more common 
than previously thought, especially in eukaryotic cells;^ for instance, it has been estimated 
that as much as 17% of the proteins in Drosophila are wholly unstructured.— In many cases, 
intrinsically unstructured proteins adopt specific structures upon binding to their biological 
targets.^* 3 - Folding and binding are then coupled, thus establishing a direct link between 
folding and function. 

It has been suggested that unstructured proteins offer several advantages in cellular 
regulation.^ Being unstructured might, for instance, allow one protein to interact with 
several targets. It might also be useful for control purposes due to rapid turnover. Further- 
more, it has been argued that being unstructured might facilitate the binding of the protein 
to a target,— by increasing the 'capture radius'. This mechanism, termed 'fly casting', 5 
was analyzed using a low- dimensional representation of the binding process. More recently, 
the thermodynamics of coupled folding-binding were examined using a generalized random 
energy model.— 

Here we study coupled folding-binding by computer simulations of a simple chain-based 
model. Simulating coupled folding-binding of a chain in a controlled manner requires proper 
sampling of the full conformational space and is harder than simulating folding of an isolated 
chain, since rigid-body translations and rotations must be taken into account. The complex- 
ity of the problem makes it highly desirable to study coarse-grained models before entering 
high-resolution modeling. For the present study, we use the minimal two-dimensional hy- 
drophobic/polar HP lattice model of Lau and Dill.— This model has been widely used to 
investigate basics of protein folding. 8 This and similar lattice-based models have also been 
used to gain insights into topics such as protein evolution,— i 10 ' 11 prion-like conformational 
propagatio n 1 ^ 13 and protein aggregation.^*^ 

Our study is inspired by recent experiments by Wahlberg et alM on the in vitro evolved 
protein Zspa-i- This sequence was engineered^ from the Z domain of staphylococcal pro- 
tein A, a well characterized three-helix-bundle protein.— It was selected for binding to the 
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Z domain itself. In the Z:Zspa-i complex, Zspa-i adopts a structure similar to the solution 
structure of the Z domain*!* However, iu solution, Z SM _, turns out to be structurally 
disordered.— The engineered ZgpA-i thus exhibits coupled folding-binding. In Ref. |2C| com- 
puter simulations of the solution behaviors of ZgpA-i and the Z domain were performed, 
using a relatively detailed off-lattice model with 5-6 atoms per amino acid.— Simulating the 
binding behavior of an unstructured protein like ZgpA-i at this level of resolution remains, 
however, a challenge. 

The present study consists of two parts. First, we study binding statistics for very short 
HP chains with up to iV = 14 monomers, to get an idea of how likely binding is to occur for 
stable and unstable sequences, respectively, in this model. We then perform a more detailed 
study of two iV = 25 sequences. In particular, this study allows us to compare the free 
energy of coupled folding-binding with free-energy profiles for docking and for folding of an 
isolated chain. 

II. MODEL AND METHODS 

In the HP model,— the protein chain is represented by a string of hydrophobic (H) or 
polar (P) beads on the square lattice. Adjacent beads along the chain are connected by 
links of unit length. It is forbidden for two beads to occupy the same lattice site. Two 
beads that are neighbors on the lattice but non- adjacent along the chain are said to be in 
contact. The energy that a configuration gets is determined by the number of HH contacts, 
each HH contact being assigned an energy e < 0. This defines the model for a single chain 
in isolation. 

Here we study single chains interacting with some fixed target structure. Two types of 
targets are considered. In Sec. 1111 Al the target is an immobilized HP chain, and the inter- 
chain interactions are taken to be the same as the intra-chain interactions; that is, each 
HH contact between the target and the flexible chain is assigned the energy e. Any cross 
HH contact is given this energy whether or not it is present in the final bound state. In 
Sec. IIIIB1 the target is an extended structure with one particular binding site. In this case, 
there is only one specific bead of the flexible chain that can interact with the binding site of 
the target. A contact between these two sites is taken to be favorable by an energy e\ } = 3e. 

To study the thermodynamics of these systems, we use Monte Carlo methods. The 
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moving chain is allowed four types of moves: local one- and two-bead moves, non-local pivot 
moves, and one-step translation moves (left, right, up or down). All moves are subject to 
a Metropolis accept/reject step, to ensure that detailed balance is fulfilled. One round of 
simulation, called one 'sweep', consists of (at most) N — 1 one-bead moves, N — 2 two-bead 
moves, one pivot move and one translation move, N being the number of beads. 

In Sec. IIII B| in addition to the thermodynamic simulations, we also study the binding 
process as a function of Monte Carlo time. These calculations follow exactly the same proto- 
col, except that the non-local pivot update is omitted, to avoid large unphysical deformations 
of the chain. 

With its simplified conformational space and its minimal two-letter alphabet, the HP 
model is not meant for studies of specific proteins, but rather to shed light on general ques- 
tions for generic sequences. It is worth noting 8 that if the principle of minimal frustration 21 
holds, then the native structures of functional proteins should be strongly favored by the 
hydrophobicity pattern alone. This suggests that the HP model, despite its simplicity, might 
be able to capture non-trivial features of the mapping from sequence to native structure. 
A statistical analysis of HP model sequences with unique ground states lends support to 
this view; it turns out that the hydrophobic beads are anticorrelated along the chains , 22123 
which is the same behavior that real (globular) protein sequences show . 22 1 24 A recent study 
of the distribution of hydrophobicity in protein sequences can be found in Ref. |25[ For a 
study of designed hydrophobic/polar copolymers with positive hydrophobicity correlations, 
see Refs. QH 

The HP model has the advantage over high-resolution models that exact results can be 
obtained for short chains, by exhaustive enumeration of all possible states. On the square 
lattice, all possible sequences with unique ground states (in isolation) have been determined 
for iV < 25,— along with the corresponding structures. Below, we make use of these results. 
Sequences having a unique ground state will be referred to as stable. 
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III. RESULTS AND DISCUSSION 



A. Statistics for Short Chains 

To get an idea of how likely different binding behaviors are to occur in this model, we 
start by studying a large number of short sequences with up to iV = 14 beads. Each 
sequence interacts with some immobilized target chain with the same length. All target 
sequences considered are stable and held fixed in their unique ground-state conformations. 
The two-chain system is contained in a box of size 2N x 2N. The target is fixed in the 
middle of the box so that the moving chain can attack it from any side. We say that the 
moving protein binds to the target if there is a unique minimum-energy configuration for 
the two-chain system. Whether or not this criterion is met can be determined with high 
confidence by Monte Carlo methods for these chain lengths. As an operational criterion for 
binding we require that the minimum-energy configuration is unique and is visited at least 
10 times during the course of the simulation. To ensure that these visits to the bound state 
are 'independent', a visit is counted only when the system comes to this state after going to 
some state which had at least 20% higher energy than the bound state. 

We first test all possible sequences with a given N for binding to one particular target 
conformation with that N, using N = 10, 11 and 12. As target, we somewhat arbitrarily 
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pick the first entry in the list of all stable ground states obtained in Ref. 
of these calculations are summarized in Table HJ From this table it can be seen that the 
fraction of binding sequences is a few per cent, which is comparable to the number of stable 
sequences. A clear majority of the binding sequences are unstable, which means that docking 
cannot be the mode of binding. That most binding sequences are unstable does not mean 
that such sequences have an intrinsically higher propensity to bind than stable ones, but 
rather it merely reflects the fact that most sequences are unstable. 

As mentioned in the introduction, one suggested use of being unstructured is that it might 
enable the chain to bind to different targets. In the model, there are sequences showing this 
behavior. An example of this can be found in Fig. ^ which shows an unstructured N = 12 
sequence that is able to bind to two different targets (both of which are held fixed). The 
structure that the unstable sequence adopts upon binding is seen to depend on the target. 

Finally, restricting ourselves to stable sequences, we test for self-binding, by using the 
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FIG. 1: Same sequence (HHPHHPPHPPHP) binding to two different targets. The moving se- 
quence (thin lines) adopts different folds as it binds to the two fixed targets (thick lines). 

same sequence as target and moving chain simultaneously. In obtaining Table H] above, this 
calculation was done for those three sequences that served as targets. Using N = 12 and 
N = 14, we now consider all stable sequences. Table ITT1 shows the results of these calculations. 
We see that ~10% of the stable sequences are self-binding both for N = 12 and N = 14. For 
N = 12, we observe that all the 9 self-binding stable sequences bind in their ground-state 
conformations. For N = 14, on the other hand, there exist 2 self-binding stable sequences 
for which the bound structure differs from the isolated one, while the remaining 40 keep the 
same structure after binding. 

Let us briefly summarize the results presented so far. Our study of all possible sequences 
for N = 10, 11 and 12 shows that a significant fraction (a few per cent) can bind to 
the targets considered. Furthermore, most of these binding sequences are unstructured in 
isolation. Our study of all stable sequences for N = 12 and N = 14 shows that ~10% of 



TABLE I: Numbers of sequences for iV = 10, 11 and 12 that can bind to the ground-state 
conformations of HHPPHPPHPH, HHPPHPPPPHP and HHPPHPPHPHPH, respectively. The 
stable sequences for these N are known from previous work i 23 i 29 



N = 10 



JV = 11 



N = 12 



Total no. of sequences, 2^ 
Stable sequences 3. 
Binding stable sequences' 3 
Binding unstable sequences 
Total no. of binding sequences 



1024 
6 

20 
20 



2048 
62 
6 
44 
50 



4096 
87 
18 
140 
158 



a A sequence is stable if it has a unique ground state in isolation. 

b A sequence is binding if the system of moving chain plus target has a unique ground state. 
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them are self-binding. For most of these self-binding sequences, the complexed conformation 
is identical to the stable structure in isolation, indicating the possibility of a docking-like 
binding behavior. 

B. Coupled Folding-Binding versus Docking 

Having studied binding propensities of short chains, we now turn to a more detailed study 
of two longer chains with iV = 25, whose binding behaviors are contrasted. This time the 
target is an extended structure located in the left-bottom end of an iV x iV box, in which the 
chain is confined (see Figs. [21 and EJ). The target has only one binding site, and the binding 
is specific in that only one bead of the moving chain can interact with this binding site. 

The two sequences considered are given in Table ITTT1 The first sequence, called S, was stud- 
ied in isolation in Ref. |3- It was obtained by applying a sequence optimization algorithm 3 ^ 
to 326 stable sequences all having the same unique ground-state conformation. It turns out 
that S has a ground-state energy of 13e. The next most favorable conformations have only 
11 HH contacts, so there is an energy gap of 2e, making S unusually stable. The other 
sequence studied, called U, is unstable. Its degenerate ground state has an energy of lie. In 
our calculations, we observed five different conformations of U with this energy. Sequence U 
is a close analog of sequence S with just one hydrophobic bead mutated to a polar one. The 
two sequences were taken to be similar in order for their bound structure to be the same. 

Our thermodynamic simulations of these two sequences were started from random con- 
figurations and contained 10 8 Monte Carlo sweeps each. The results show first of all that 

TABLE II: Numbers of stable sequences that can bind to themselves for N = 12 and N = 14. The 
total numbers of stable sequences are from previous work. 23 ' 29 The last row gives the numbers of 



self-binding stable sequences that bind in their ground-state conformations. 

N = 12 N = 14 

Total no. of sequences 4096 16384 

Stable sequences 87 386 

Self-binding stable sequences 9 42 

Self-binding stable sequences with unchanged structure 9 40 



7 



FIG. 2: Typical binding behavior of sequence S. The chain first folds and then just translates to 
bind to the target. Filled and open circles represent hydrophobic and polar beads, respectively, 
and the circle with a cross represents the binding site. 




FIG. 3: Snapshots from a simulation of sequence U. After binding to the target, the chain rear- 
ranges itself into the minimum-energy state (right-most image). Symbols are as in Fig. [3 

both these sequences do bind to the chosen target, in the sense that the system of target 
plus moving chain has a unique minimum-energy state. This state, referred to as the bound 
state, was visited many independent times in the simulations. The bound structure is the 
same for both chains (see Figs. El and El), and coincides with the stable structure of S in 
isolation. The energy of the bound state is 13e + eb for S and lie + eb for U. 

To characterize the binding behaviors of these sequences, we monitor the following two 

TABLE III: The two N = 25 sequences studied, S and U. The position at which they differ is 
underlined. The bead that can interact with the target is in bold type. 

S: HHHHH PPHPP HPHPH PPHPH PHPHP 
U: HHHHH PPHPP HPHPH PPPPH PHPHP 



S 



quantities: 



1. The binding parameter /. / is the geometric distance between the binding site of the 
target and that bead of the moving chain that can interact with the target. As binding 
progresses, the value of I reduces. 

2. The folding parameter Q' . The chains studied form 14 internal contacts, called native 
contacts, in their bound states. The number of these native contacts being present, 
Q, provides a measure of the 'nativeness' of the chain. As the chain folds to its 
bound-state structure, the value of Q increases. It turns out that the value Q = 13 is 
impossible to attain for these sequences. In our free-energy calculations, we therefore 
use a folding parameter Q' defined by 



Figures 0] and El show the free energy calculated as a function of these two vari- 
ables, F(I,Q'), for the sequences S and U, respectively. The free energy is defined by 
P(I,Q f ) = exp[—F(I,Q')/kT], where P(I,Q') is the joint I,Q' probability distribution. 
The temperature is taken as T = e/2.6k for S and T = e/3.2k for U (k is Boltzmann's 
constant). These temperatures are chosen such that the binding probabilities for the two 
chains are significant and close to each other (« 35%). Note that the systems are in their 
minimum-energy states if and only if (J, Q') = (0, 13). 

For sequence S (Fig. HJ), we find that F(I, Q') has a simple shape with two narrow valleys 
along the lines Q' = 13 and 1 = 0, respectively, and a broad and shallow minimum centered 
at (I,Q') ~ (11,5), where the chain is unbound and unfolded. This suggests that there are 
two very different major binding modes for this sequence. One way for the chain to reach 
its bound state is along the line Q 1 = 13, which corresponds to rigid docking; the chain first 
folds and then moves towards the binding site. The other major binding mode is along the 
1 = valley. Here the chain first attaches to the binding site and then folds to its final 
shape. Following Ref. 0, we refer to this behavior as the fly-casting mechanism. 

Let Pd and Pf denote the probabilities of finding the system in the Q' = 13 and 1 = 
corridors, respectively (not counting the bound state); that is, 
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FIG. 4: Free energy F(I, Q') for sequence S at T = e/2.6k. The bound state, in which the chain is 
both folded and attached to the target, corresponds to the upper left corner of the figure. A chain 
that is folded but not attached to the target moves along the upper edge of the figure (docking), 
whereas a chain that is attached to the target but incompletely folded moves along the left-hand 
edge (fly casting). The contours are spaced at intervals of 1 kT and dark tone corresponds to low 
free energy. Contours more than 8 kT above the minimum free energy are not shown. 

where the subscripts d and f refer to docking and fly casting. For sequence S, we find that 
Pd is about twice as large as Pf (see Table ITVj) . suggesting that docking is the preferred 
binding mode for this sequence. 

For the unstructured sequence U (Fig. [HJ), there is no free-energy valley corresponding 
to rigid docking, so Pd is small (see Table I1VJ) . There is, by contrast, an I = valley 
corresponding to fly casting for this sequence, too. The population Pf is, in fact, higher 
for U than for S (see Table HV)) . The 1 = valley is, for both sequences, separated from 
rest of the conformational space by a free-energy ridge. For sequence U, there are some 



TABLE IV: The probabilities Pd an d Pf as defined by Eq. © a along with the binding probability 
Pb = P(I = 0, Q' = 13), for the sequences S and U. 

S U 

P b 0.35 0.35 

P d 0.32 0.05 

P f 0.17 0.23 

a The subscripts d and f refer to docking and fly casting. 
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FIG. 5: Same as Fig. 0]for sequence U at T = e/3.2k. 

narrow troughs in these hills towards the 1 = corridor, which should make it easier for 
this sequence to reach this corridor. In addition to the 1 = valley, sequence U exhibits 
two somewhat less pronounced valleys along the lines Q' = 9 and Q' = 12, respectively, as 
well as a broad unbound and unfolded minimum at lower Q'. Note that there exist Q' = 12 
conformations with minimal intra-chain energy for this sequence. The presence of the Q' = 9 
and Q' = 12 valleys suggests that this chain can follow several different paths to its bound 
state. If the chain follows one of the fixed-Q' valleys, it reaches the binding site without 
having its full bound-state structure; folding is completed after the chain has attached itself 
to the target. When increasing the box size to 3iV x 3iV, we observed valleys at the same 
values of Q' which were stretched along I. 

It is instructive to take a closer look at the free energy in the docking and fly-casting 
corridors, respectively. Figured shows the free energy along the docking corridor, F(I, Q' = 
13), for sequence S. We see that in order to reach the bound state, the chain has to pass a 
free-energy barrier with a height of ~ 3 kT. This barrier is entropic. Figure \§jp shows the 
free energy along the fly-casting corridor, F(I = 0,Q'), for both sequences. Two features 
are worth noting. First, the shape of the curve is roughly the same for the two sequences, 
although the curve is shifted downwards for sequence U which spends more time in this 
corridor. Second, there is no major free-energy barrier in this corridor (the highest barrier 
is about lkT), and the last part towards the bound state (from Q' = 8-9 to Q' = 13) is 
downhill in free energy. 

It is also interesting to compare the free energy in the fly-casting corridor (Fig. EJd) with 
the free energy of folding for sequence S in isolation. For this purpose, we performed a sim- 
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FIG. 6: Free-energy profiles for the docking and fly-casting corridors, (a) F(I, Q' = 13) for 
sequence S. (b) F(I = 0,Q') for sequence S (dotted line) and sequence U (solid line). The arrows 
indicate the location of the bound state, in which the chains are both folded and attached to the 
target. The temperature is e/2.6k for S and e/3.2k for U (same as in Figs. |1] and 03 respectively). 
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FIG. 7: Two free-energy profiles for sequence S. The dashed line shows the free energy along 
the fly-casting corridor, F(I = 0,Q'), at T = e/2.6k (same as in Fig. H>p)- The full line shows 
F(Q') for the isolated chain at T = e/2.268/c. U and F refer to unfolded and folded conformations, 
respectively. 

ulation of S in isolation at T = e/2.268/c, where the native population (29%) is comparable 
to the binding probability at T = e/2.6k (35%). From Fig. it can be seen that the free 
energy for the isolated chain is markedly different from that for the fly-casting corridor. In 
particular, we see that the bound-state minimum of the interacting chain is broader than 
the native minimum of the isolated chain; the isolated chain must reach Q' = 11 before the 
free energy starts to decrease. 
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We also performed Monte Carlo-based kinetic simulations (see Sec. |TTJ) of the sequences 
S and U, using the same temperatures as in the thermodynamic runs. For each sequence, 
we carried out a set of 500 simulations, each containing 5 x 10 6 Monte Carlo sweeps. The 
simulations were started from configurations obtained by short runs at a higher temperature 
(T = e/2k). The resulting configurations from these preparatory runs (the starting config- 
urations for the kinetic runs) showed a wide variation in both location and conformation 
for the moving chain. The averages of / and Q over this ensemble were about 12 and 4-5, 
respectively. 

Figure |H1 shows the evolution of this ensemble of 500 systems with Monte Carlo time in 
the (J, Q) plane (not Q'). This parametric plot, with Monte Carlo time as a parameter, has 
the advantage of being independent of the difference in time scale between the two different 
temperatures used. To reduce noise, each data point represents an average over 10 4 Monte 
Carlo sweeps, thus giving 500 data points in total. In Fig. we also indicate the equilibrium 
values of / and Q, which were obtained by a separate, very long simulation. The equilibrium 
values are given by ((/), (Q)) = (6.1, 12.1) for S and ((I), (Q)) = (4.9, 11.1) for U. Note that, 
as expected, S has a higher value of (Q) and U has a lower value of (/). 

For sequence S, we see that the relaxation towards the equilibrium point has a clear two- 
step character, where the first step corresponds to folding (increasing Q) and the second 
step to binding (decreasing /). This behavior is indeed what one expects if docking is the 
preferred binding mode. For sequence U, there is no such clear distinction between folding 
and binding; I and Q evolve in a more correlated manner. Hence, the results of the kinetic 
simulations support the conclusions from the free-energy analysis. 

IV. SUMMARY 

We have studied different binding mechanisms of polypeptide chains using the minimal 
HP model. The aim of this study is to understand the nature of an important aspect of 
protein-protein recognition and it can serve as a guideline for further, detailed studies of 
this complex phenomenon. 

For small chain lengths, we found that the fraction of sequences that can bind to a given 
(fixed) target is comparable to the fraction of sequences that are stable in isolation. The 
overlap between these two sets of sequences is small; most binding sequences are unstable 
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FIG. 8: Parametric plot of Q against / with Monte Carlo time as a parameter for sequence 
S (dotted line) and sequence U (solid line). The curves represent averages over 500 simulations 
and the bottom-right ends correspond to time zero. The unfilled and filled squares represent the 
equilibrium points ((Q),(I)) for S and U, respectively. Temperatures are as in Figs. |I] and [SJ 

in isolation and fold upon binding to the target. 

We then compared the binding behaviors of two related N = 25 sequences. One of the 
sequences, S, had been optimized for high stability in isolation, 23 and was found to have 
rigid docking as its preferred binding mode. The other, slightly mutated sequence, U, was 
found to have a bound structure identical to that of S, but to be unstructured in isolation. 
In these respects, sequence U is reminiscent of the engineered protein Z SP a_i mentioned in 
the introduction. It turned out that sequence U tends to first attach to the target and then 
fold, which supports the conclusion^ that unstructured chains prefer binding through a fly- 
casting mechanism. The free-energy profile associated with this binding mode was found to 
lack high barriers and to exhibit a broad bound-state minimum, compared with the native 
minimum of sequence S in isolation. The observed difference in shape between these free 
energies suggests that by attaching itself to the target, the chain becomes able to fold more 
efficiently than it does in isolation. The free energy of docking obtained for sequence S 
contains, as expected, a significant entropic barrier, which the chain has to overcome in 
order to reach the bound state. 

These free-energy profiles were obtained using two particular sequences in this model. 
Therefore, it should be stressed that the shape of the free energy was not considered at 
all when selecting these two sequences and the target structure; the goal was just to have 
two binding sequences with the same bound structure but different stability properties in 
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isolation. As a result, we believe that the trends seen have some generality. Nevertheless, 
it is clear that it would be highly desirable to extend these calculations to other and longer 
sequences and to more realistic off-lattice models. Furthermore, it would be interesting to 
study the kinetics over a larger time interval; this requires a larger number of systems than 
500 which we used, in order to keep the statistical errors under control. 
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